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Summary 

In  this  report  we  summarize  numerical  results  of  supersonic  flows  over  open  cavities  in 
the  setting  of  a  dual-mode  ramjet/scramjet  engine.  To  calculate  the  unsteady  cavity  flows,  we 
employ  the  Space-Time  Conservation  Element  and  Solution  Element  (CESE)  method,  a  novel 
numerical  method  based  on  a  unified  treatment  of  space  and  time  in  calculating  flux  balance. 
Supersonic  cavity  flows  with  and  without  fuel  injection  are  studied  to  understand  the 
mechanisms  of  mixing  enhancement  and  flame  holding  by  cavities.  Without  injection, 
numerical  results  compared  favorably  with  the  experimental  data  for  dominant  frequencies  and 
time-averaged  pressure  coefficients  inside  the  cavities.  With  an  upstream  injection,  the  flow 
oscillations  are  drastically  suppressed.  In  a  downstream  injection  arrangement, 
cavity-generated  acoustic  waves  and  vortices  greatly  enhance  fuel/air  mixing.  Numerical 
results  show  that  the  CESE  method  provides  high-fidelity  numerical  results  of  unsteady  flows 
in  the  advanced  scramjet  engine  concept. 

The  results  reported  here  have  been  published  in  the  following  journal  and  conference 


papers: 

1 .  C.-K.  Kim,  S.-T.  J.  Yu,  and  Z.  Zhang,  “Direct  Calculation  of  High-Speed  Cavity  Flows 
in  a  Scramjet  Engine  by  the  CESE  method,”  AIAA  J.,  vol.  42,  no.  5  (2004)  pp.  912-919. 

2.  Z.  Zhang,  S.-T.  J.  Yu  and  S.-C.  Chang,  “A  Space-Time  Conservation  Element  and 
Solution  Element  Method  for  Solving  the  Two-  and  Three-Dimensional  Euler  Equations 
by  Quadrilateral  and  Hexahedral  Meshes,”  Journal  of  Computational  Physics,  vol.  175, 
no.  1  (2002)  pp  168-199. 

3.  C.-K.  Kim,  K.-S.  Im  and  S.-T.  J.  Yu,  “Numerical  Simulation  of  Unsteady  Flows  in  a 
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Scramjet  Engine  by  the  CESE  Method,”  the  JANNAF  Meeting,  Aero  Propulsion,  June 
2002,  Destin,  FL. 
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1.  Introduction 


Fuel  injection,  ignition,  and  flameholding  are  challenging  issues  for  high-speed  combustion.  In 
a  viable  scramjet  engine,  the  fuel  injection  method  employed  must  provide  rapid  fuel/air 
mixing  with  minimum  total-pressure  loss  in  the  airstream.  A  stable  flameholding  system  under 
a  wide  range  of  operating  conditions  is  critical  to  sustain  the  supersonic  combustion.  Recently, 
cavity-based  flameholders,  an  integrated  mixing-enhancement  and  flameholding  approach, 
have  attracted  considerable  attention  in  the  scramjet  community.  Under  suitable  conditions, 
flow  recirculation,  or  the  trapped  vortices,  significantly  increases  the  flow  residence  time  of  the 
fluid  entering  the  cavity.  A  pilot  flame  could  be  set  up  inside  the  cavity  to  provide  a  pool  of  hot 
chemical  radicals,  which  in  turn  would  reduce  the  ignition  delay  of  the  air/fuel  mixture  in 
airstream  and  thus  sustain  high-speed  combustion. 

High-speed  cavity  flows  are  inherently  unsteady,  involving  both  broadband  small  scale 
fluctuations  typical  of  turbulent  flows,  as  well  as  distinct  resonance  with  harmonic  properties  in 
its  frequencies  and  amplitudes.  In  the  past,  it  has  been  demonstrated  that  the  aspect  ratio  of  the 
cavity  and  free  stream  flow  conditions  are  the  critical  parameters  dominating  the  complex  flow 
features,  including  boundary  layer  separation,  compressible  free  shear  layer  with  shedding 
vortices,  linear/nonlinear  acoustic  waves,  and  complex  shock  and  expansion  waves  interacting 
with  vortices  and  acoustic  waves. 

In  the  setting  of  wheel  wells  and  bomb  bays,  previous  studies  for  high-speed  cavity 
flows  showed  that  cavity  flows  could  be  categorized  into  the  following  two  groups:  (i)  open 
cavity  flows,  when  L/D  <  7~10,  and  (ii)  closed  cavity  flows,  when  L/D  >  7—10,  where  L 
denotes  the  length  of  the  cavity,  and  D  the  depth.  In  flows  over  cavities  of  large  aspect  ratios 
(L/D  >  7—10),  the  separated  free  shear  layer  emanating  from  the  upstream  corner  of  the  cavity 
reattaches  to  the  bottom  wall  of  the  cavity  and  results  in  two  separated  recirculation  zones  near 
the  two  comers  between  the  lateral  walls  and  the  cavity  floor.  The  resultant  low  pressure  zones 
at  the  lower  comers  and  high  pressures  on  the  cavity  floor,  where  the  shear  layer  reattaches, 
lead  to  significant  drag  and  pressure  loss  of  the  airstream.  In  this  case,  mass  addition/ejection 
into/from  the  cavity  by  aerodynamic  unsteadiness  is  low  to  moderate,  and  the  flow  is  referred  to 
as  “closed.” 
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On  the  other  hand,  flows  over  cavities  with  smaller  aspect  ratios,  L/D  <7-10,  result  in 
reattachment  of  the  free  shear  layer  to  the  rear  bulkhead  of  the  cavity.  The  impingement  of  the 
free  shear  layer  on  the  rear  lateral  wall  causes  violent  unsteady  motions  and  results  in 
significant  periodical  mass  addition/ejection  near  the  rear  bulkhead  of  the  cavity,  and  the  flows 
are  referred  to  as  “open.”  The  wave  patterns  of  open  cavity  flows  could  be  further  categorized 
into  (i)  transverse  mode  for  very  short  cavities,  L/D  =  1,  and  (ii)  longitudinal  mode  for  longer 
cavities,  e.g.,  2-3  <  L/D  <  7-10.  In  short  cavities,  L/D  <  2,  only  one  main  vortex  inside  the 
cavity  is  sustained  by  the  driving  shear  layer  spanning  the  top  of  the  cavity.  The  up  and  down 
motions  of  the  single  main  recirculation  bubble  generate  acoustic  waves,  which  by  and  large 
propagate  in  the  direction  perpendicular  to  the  free  shear  layer,  provided  the  free  stream  is 
transonic.  The  propagating  waves  are  referred  to  as  in  a  transverse  mode.  On  the  other  hand, 
when  the  cavity  is  longer,  2-3  <  L/D  <7-10,  multiple  moving  vortices  occur  inside  the  cavity 
leading  to  complex  interactions  among  trapped  vortices,  propagating  and  rebounding  pressure 
waves,  and  the  flapping  free  shear  layer.  In  general,  the  rebounding  pressure  waves,  while 
interacting  with  the  free  shear  layer,  drastically  amplify  the  growth  rate  of  the  free  shear  layer, 
which  in  turn  sheds  enormous  vortices  propagating  towards  and  impinging  on  the  aft  wall  of  the 
cavity.  Due  to  propagating  vortices  in  the  streamwise  direction  and  the  rebounding  pressure 
waves,  prevalent  acoustic  waves  propagate  in  the  longitudinal  direction  outside  the  cavity  into 
the  downstream  area.  If  the  airstream  is  transonic  or  subsonic,  the  acoustics  would  transversely 
propagate  into  the  upstream  areas. 

In  the  setting  of  supersonic  combustion  inside  a  scramjet  engine,  trapped  vortices  inside 
cavities  could  be  useful  for  flame  holding.  Moreover,  cavity  resonance,  which  produces 
periodic  mass  addition/expulsion  with  large  flow  structures,  could  be  useful  for  mixing 
enhancement.  Simultaneously,  cavity  drag  must  be  minimum,  e.g.,  much  less  than  that  of  a 
bluff  body,  and  thus  only  causes  acceptable  pressure  loss.  Gruber  et  al.  [1, 2]  have  developed  a 
dual-mode  ramjet/scramjet  engine  concept,  which  is  envisioned  to  use  hydrocarbon  fuels  for  a 
flight  regime  of  Mach  numbers  from  3  to  6  -  9.  In  their  supersonic  combustors  [1,  2],  open 
cavities  with  aspect  ratios  about  5  <  L/D  <  8  have  been  tested  in  conjunction  with  various  fuel 
injection  schemes.  Numerical  simulation  of  cavity  flows  has  been  conducted  by  Baurle  et  al. 
[3].  The  results  showed  that  the  cavities  have  great  potential  to  be  a  viable  combined 
flameholder/mixing  enhancement  device  for  scramjet  engine  combustor.  Similar  ideas  have 
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also  been  independently  proposed  and  tested  by  Yu  et  al.  [4],  In  particular,  Yu  et  al.  [4]  has 
tested  supersonic  flows  passing  multiple  cavities.  Some  of  recent  results  have  been 
summarized  by  Ben-Yaker  and  Hanson  [5]. 

In  the  past,  extensive  experimental  and  theoretical  studies  on  cavity  flows  have  been 
conducted  for  applications  in  wheel  well  and  bomb  bay,  and  flow  characteristics  such  as  the 
oscillation  frequency  and  amplitudes  at  various  locations  in  the  cavity  have  been  reported 
[6-11].  However,  it  is  difficult  to  directly  apply  this  knowledge  base  to  cavity  flows  for  the 
advanced  scramjet  engines  due  to  much  shorter  length/time  scales  in  scramjet  engine. 
Additional  complexity  associated  with  fuel  injection  also  warrants  further  studies  because 
cavity  flows  and  the  associated  acoustics  would  be  drastically  changed  by  the  fuel  injection 
schemes  employed.  In  particular,  inherent  oscillations  of  cavity  flows  may  be  significantly 
suppressed  by  an  upstream  injection  [4,  9-11]. 

In  the  present  report,  we  will  focus  on  time-accurate  calculation  of  supersonic  cavity 
flows  in  the  setting  of  a  dual  mode  ramjet/scramjet  engine  combustor  [1,  2].  The  objectives  of 
the  present  study  are:  (i)  to  validate  the  numerical  results  by  assessing  the  calculated 
frequencies  and  amplitudes  of  pressure  oscillations  and  compared  them  with  previous  reported 
data;  (ii)  to  assess  the  fuel/air  mixing  enhancement  based  on  applying  upstream  as  well  as 
downstream  injection  to  cavity  flows;  and  (iii)  to  demonstrate  the  capabilities  of  the  CESE 
method  for  capturing  complex  flow  features  of  the  supersonic  cavity  flows. 

The  rest  of  the  report  will  be  organized  as  follows.  Chapter  2  reviews  the  model 
equations  to  be  solved  by  the  CESE  method.  Detailed  derivation  for  the  flow  equations  of  the 
gas  mixtures  is  presented.  Chapter  3  provides  background  information  of  the  CESE  method. 
Chapter  4  shows  numerical  solutions,  including  comparison  between  the  numerical  results  and 
previously  reported  data.  Moreover,  we  will  show  the  effects  by  both  upstream  and  downstream 
injection  on  pressure  oscillations,  acoustics,  and  vortices,  leading  to  effects  on  fuel/air  mixing 
and  flameholding.  We  then  offer  concluding  remarks  and  provide  cited  references. 
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2.  Governing  Equations 


In  this  chapter,  comprehensive  discussions  are  provided  for  the  governing  equations  of  gas 
mixtures.  Although,  two-dimensional  solutions  will  be  provided  in  Chapter  3,  we  will  focus  on 
one-dimensional  equations  with  detailed  derivation.  Extension  to  multi-dimensional  flow 
equations  is  straightforward  and  will  only  be  briefly  illustrated  at  the  end  of  this  chapter. 

2.1  Mixtures  of  Thermally  Perfect  Gases 

Consider  a  mixture  of  N  different  thermally  perfect  gas  species,  which  generally  are  not 
calorically  perfect,  i.e.,  the  specific  heats  of  each  species  is  not  constant.  Let  the  species  be  in 
thermal  equilibrium  with  each  other.  As  such  all  species  share  the  common  mixture  absolute 
temperature  T and  each  individual  gas  species  /  (/  =  1,2,..., TV)  satisfies  the  perfect  gas  law 

P,=PAT  (2.1) 

Here/’ ,  pj ,  and  R  are  the  partial  pressure,  the  mass  density,  and  the  gas  constant  of  species  i, 
respectively.  Note  that Rj=RJ1/Mj,i  =  1,2,...,  N ,  with  Ru  and  M.,  respectively,  being  the 

universal  gas  constant  and  the  molecular  weight  of  species  /.  In  addition,  it  is  assumed  that  the 
static  pressure  p  of  the  mixture  can  be  determined  using  Dalton’s  law,  i.e., 

P,=t,P,=Tfjp,R,  (2.2) 

i=l  i=l 

To  proceed,  a  set  of  definitions  is  given  in  the  following.  The  mass  fraction  of  species  i  is 

p 

Yt  =  (2.3) 

P 

where  p  is  the  mass  density  of  the  mixture,  i.e., 

def  N 

P  =  Zp,  (2-4) 

7=1 

Note  that,  by  definition, 

(2.5) 

i=\ 

For  a  reason  to  be  shown  (see  Eq.  (2.14)), 
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(2.6) 


def  N 
i= 1 

is  referred  to  as  the  gas  constant  of  the  mixture.  Let  ei  and  s: ,  respectively,  denote  the  specific 

internal  energy  and  the  specific  entropy  of  species  i.  Then  the  specific  internal  energy  and  the 
specific  entropy  of  the  mixture  are 

def  N 


and 


=5>, 


def  N 
1=1 


(2.7) 


(2.8) 


respectively.  Note  that,  because  the  specific  internal  energy  of  any  thermally  perfect  gas 

*/  jp  / 

species  is  a  function  of  temperature  only,  each  et  is  a  function  of  T  only.  Let  e,  =  for 

each  Then 


def  N  • 

Cv  =  ^fie'  >0 

/=1 


(2.9) 


Replacing  YN  by  ^  ,=i  %  (See  Eq.  (2.5)),  Eq.  (2.9)  implies  that 

.  N- 1  .  . 

cv  =eN  +  IlXft  — sn)  >  0 


(2.10) 


/=i 


Using  Eqs.  (2.3)  and  (2.4),  and  replacing  pN  by  pt  (See  Eq.  (2.4)),  Eq.  (2.2),  (2.6) 

and  (2.7)  imply  that 


N- 1 


P=  pRn+T.Pi(Ri-Rio  T 

/=! 

(2.11) 

,=l  P 

(2.12) 

N~l  p 

e  ~  eN  "*■  2]  ( ei  ~  eN  ) 

,=I  P 

(2.13) 

respectively.  The  key  conclusions  that  can  be  drawn  from  Eqs.  (2.1 1)-(2.13)  are  given  in  the 
following.  Eqs.  (2.1 1)  and  (2.12)  imply  that 
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p  =  pRT 


(2.14) 


i.e.,  R  is  the  gas  constant  of  the  mixture.  Note  that  Ri ,  i  =  1, 2  are  assumed  to  be  constant 
in  the  current  study.  As  such,  with  the  aid  of  Eq.  (2.12),  Eq.  (2.1 1)  implies  that 

dp  =  pRdT  +  TRNdp  +  TY  (R,  ~  RN)dPi  (2.15) 

i= 1 

Because  each  ei  is  a  function  of  T  only,  Eq.  (2.13)  implies  that  e  is  a  function  of  T,  p  and  p, , 
i  -  1,2,...,  jV-1  .  As  such,  with  the  aid  of  Eq.  (2.10)  and  the  relation  de:  =  etdT ,  one  has 


de  =  c,dT-\!Zp,(ei-en)\^+-fjiel-e„)dp,  (2.16) 

M  P  Pi=\ 


It  follows  from  Eq.  (2.16)  that 


(2.17) 


P’Pl.Pl  y--yPN-\ 


Because  holding  the  value  of  pand  pn  i  =  1,2,...,  A -1  constant  is  equivalent  to  holding  the 
value  of p, ,  i  =  1,2,..., N  constant  (see  Eq.  (2.4)),  Eq.  (2.17)  implies  that  cv  is  the  specific 


heat  of  the  gas  mixture  at  constant  volume.  Because  e  is  a  function  of  T, 
p  and  pi ,  i  =  1,2,..., iV  — 1 ,  implicitly  T  is  a  function  of  e,  p  and  pn  i  =  1,2 


Substituting  this  functional  relation  into  Eq.  (2. 1 1),  one  arrives  at  the  conclusion  that  implicitly 
P  is  a  function  of  e,  p and  pt,  i  =  1,  2,  3,  ...,  iV — 1 ,  i.e.. 


p  =  pm  (e,  p,  a  ,  p2 , . . . ,  )  (2.18) 

Here  P(l)  is  a  function  of  e,  p  and p, , i  =  1,2,...,A-1. 

By  eliminating  dT  and  using  the  fact  thatcv  >  0 ,  Eqs.  (2.15)  and  (2.16)  can  be  combined  to 
yield 

#  =  — +  R„T+— dp  +  Y  (R-RN)T-~{e,-eN)  dPi  (2.19) 

L  P°v  M  J  '=1  L  Cv 

From  Eq.  (2.19),  one  arrives  at  the  conclusion  that 


r «  y  Qe’  P<P\.Pl.,P*-\ 


(2.20) 
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(2.21) 


def ,  dp .  V  T  J.  R  X’’  /•  \ 

°P  Pcv  i= 1 


And,  for  /  =  1,2,...,  A- 1, 

5p 

n  —  (— i— ) 

•“A  V  Qp'e*P*Pi,Plr.^PN-\lpi 

Here,  for  any  /  =  1,2,...,A-1, 


=  (R,_RN)T+-(e,-eH) 
c„ 


(2.22) 


'5p/ 


(¥L\ 

V  ^  '  e*P\,Plr.^PN-\  I  Pi 


denotes  the  partial  derivative  of  p  with  respect  to  pi  assuming  that,  except  p,,  all  other 
parameters  in  the  set  {e, p,p,,p2,...,pw_,}  are  held  constant.  Note  that  an  immediate  result  of 
Eqs.  (2.11),  (2.21)  and  (2.22)  is 

P^PPp+YjPiP*  (2-23) 

i=\ 

Next,  we  shall  provide  a  brief  discussion  of  the  “frozen  speed  of  sound”  of  a  gas  mixture.  As  a 
preliminary,  note  that  by  using  the  thermodynamic  relation  s,  =  st  (T,  pt ) ,  the  relation  pt-  pYt, 

and  Eq.  (2.5),  Eq.  (2.8)  implies  that  s  is  a  function  of  T,  p,  and}^  i  =  l,2,...,iV-l.  As  such 
implicitly  Tis  a  function  of  s,  p,  and  Yi,  i  =  1,2,..., N  -1 .  Substituting  this  functional  relation 
into  Eq.  (2.14)  and  observing  that  R  is  a  function  of  Yi,  i  =  1,2,...,  N only  (see  Eqs.  (2.5) 
and  (2.6)),  one  concludes  that  implicitly  p  is  a  function  of  s,  p,  and  Yi,  i  =  1, 2, ...,  N  - 1 ,  i.e., 

p  =  pm(s,p,Yl,Y2,...,YN_l)  (2.24) 


Here  P<2)  is  a  function  of  s,  p  and  Yf,  i  =  1,2,...,  N- 1 .  Let  af  be  the  frozen  speed  of  sound 


of  the  mixture,  i.e., 


(2.25) 


Note  that,  because  of  Eq.  (2.5),  holding  the  values  of  Yt  ,  /  =  1,2,...,A-1  constant 


automatically  implies  that  the  value  of  YN  is  also  held  constant.  As  such,  adding  YN  to  the  list 

of  the  subscripts  of  the  partial  derivative  in  Eq.  (2.25)  will  have  no  impact  on  the  derivative.  In 
the  following,  it  will  be  shown  that 
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a,=.  Rt{  1+— l=J-fl+— 


cj  W  P 


(2.26) 


To  proceed,  note  that  by  using  (i)  p,  =  ^ ,  (ii)  />  =  PYR,T ,  and  (iii)  the  thermodynamic  relation 


One  concludes  that 


*,=i 

T\  pf 


,  e/ 

as,  - — — 

'  T  p  Y, 


(2.27) 


With  the  aid  of  Eqs.  (2.6),  (2.9)  and  (2.28),  Eq.  (2.8)  implies  that 

A  =  £A_M£+f(Ji_w  (2.29) 

T  p  f=i 

Moreover,  with  the  aid  of  Eq.  (2.6),  Eq.  (2.14)  implies  that 

dp  =  pRdT  +  RTdp  +  pTY  R,dYt  (2.30) 

/=! 

By  eliminating  dT  and  using  the  fact  that  dYN  =  — y.._,  dYt  (see  Eq.  (2.5)),  Eqs.  (2.29)  and 
(2.30)  can  be  combined  to  yield 

dr,  (2.31) 

/>  n  .  ,  C.  C. 


(2.28) 


(2.30) 


Ml  V  C, 


An  immediate  result  of  Eq.  (2.31)  is 


=  RT(  1  +  — ) 
c„ 


(2.32) 


Eq.  (2.26)  now  follows  immediately  from  Eqs.  (2.25),  (2.32),  (2.14),  and  (2.20).  QED. 


2.2  Governing  Equations 

We  assume  that  the  N  gas  species  considered  here  share  a  common  flow  velocity  u.  Let  the  total 
specific  energy 


£  =  _ 

2 


(2.33) 
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Let  cb:,i  =  1,2  ,  be  the  net  mass  of  species  i  generated  per  unit  time  and  per  unit  volume 

from  all  chemical  reactions  involved.  Because  mass  is  conserved  during  any  chemical  reaction, 
we  have 

^&>  =  0  (2-34) 

i=\ 

Furthermore,  let 

def  def  def  def 

«I  =  p,  u2  =  pu,  u3  =  pE,  u3+i  =  Pj  (i  =  l,2,...,N)  (2.35) 

def  def  def  def 

f=pu,f=pu2  +  p,f  =  u(pE  +  p ),  ful  =  uPj  (/  =  1, 2, ... ,  N)  (2.36) 

def  def  def  def  • 

bx  =  0 ,b2=  0,  b3  =  0,  bM  =co  (/  =  1,2,..., AT)  (2.37) 

Then  the  ID  unsteady  Euler  equations  for  a  chemically  reacting  mixture  of  N  thermally  perfect 
gas  species  in  thermal  equilibrium  with  each  other  can  be  written  as 

^nL  +  ^L  =  bm  (m  =  l,2,...,N  +  2)  (2.38a) 

dt  dx 

Note  that:  (a)  In  Eq.  (2.38a),  the  equations  with  w  =  1,2,  3,  respectively,  represent  the  mass, 
momentum,  and  energy  conservation  relations  for  the  mixture.  On  the  other  hand,  those 
with  m  =  4, 5,  ...,  N  +  2 ,  respectively,  represent  the  mass  conservation  relations  for 
species /'  =  l,2,...,N-\.  Moreover,  with  the  aid  of Eqs.  (2.4)  and  (2.34),  the  mass  conservation 
relation  for  species  N  can  be  derived  from  the  mass  conservation  relations  for  the  mixture  and 
species/  =  1,2,..., /V-l .  As  such  the  former  (i.e.,  the  equation  with  m  =  N+  3)  becomes 

redundant  and  therefore  is  not  included  in  Eq.  (2.38a).  (b)  Letw ,  f ,  and  b  ,  respectively,  be  the 
column  matrices  formed  by  um,  fm,  and  bm,  m=  1, 2, ...,  N  +  2.  Then  the  matrix  form  of  Eq. 
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i.e.,  each  fm,  m=  1,2,...,  A'' +2 can  be  expressed  explicitly  in  terms  of p  and  um, 
m=  1,2,. ..,7V  +  2.  Moreover,  by  using  Eq.  (2.18)  and  the  relations 

,  p  =  ux,  pt=u3+i  (/  =  1,2,...,7V-1)  (2.40) 

which  follow  from  Eqs.  (2.33)  and  (2.35),  one  concludes  that  implicitly  p  is  also  a  function 
of  um  ,  m=  l,2,...,7V  +  2  .  Thus  implicitly  fm,  m-  1,2, 3, ...,7V +  2  ,  are  functions  of  um, 
m  =  l,2,...,7V  +  2. 

(d)  It  will  be  shown  in  Sec.  2.4  that,  implicitly®,,/  =  1,  2,  3,  ...,  TV -1  are  also  functions  of 
um,  m  —  1,2,. ..,7V  +  2. 

From  the  above  discussions,  one  concludes  that  Eq.  (2.38a)  represents  a  system  of  7V+  2 
independent  equations  for  TV  +  2  unknowns  of  um,  m  =  1, 2, ...,  TV  +  2 . 


2.3.  Jacobian  Matrix  of  Flux  Functions 


With  the  aid  of  Eqs.  (2.1 8),  (2.20)-(2.22),  and  (2.40),  an  application  of  the  chain  rule  leads  to 


$L  =  p  JlL _eW  =  dP  -_P* 

dux  p  ^  2  J  p  ’  du2  p  ’  du3  p’ 


dp 

du3+i 


(?  =  1,2,...,  TV- 1)  (2.41) 


Hereafter,  without  using  explicit  notations,  let  it  be  understood  that  a  partial  derivative  with 
respect  to  any  um  will  always  be  evaluated  assuming  that,  except  um,  all  other  independent 
variables  in  the  set  fu\,  uN+ 2}  are  held  constant. 

Let  A  be  the  (N  +2)x(N  +2)  Jacbian  matrix  with  the  element  on  the  m  th  row  and  the  ^th 


column  being  dfm  /du( .  Then  Eqs.  (2.39)  and  (2.41)  imply  that 
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0 

1 

0 

0 

0 

0  N 

Pp-— A-“2(l 
p 

2  P 

u(2-—) 

P 

Pe. 

P 

Pp, 

P* 

-  Pp,-, 

(  u  “2 

<Pp-H  +  —-pe 

2  P 

~LPe) 

P 

„  «2 

H - Pe 

P 

»0+&) 

p 

u Pp, 

uPp> 

-  UPPN-, 

u 

- Pi 

P 

RP 

P 

0 

u 

0 

0 

(2.42) 

u 

— Pi 
p 

BP 

P 

0 

0 

u 

0 

u 

- Pn-1 

Pn- i 

0 

0 

0 

u 

where 

H*E  + pi  +£■+&■)  (2.43) 

1=1  ^  Pi 

is  the  total  enthalpy  of  the  gas  mixture.  Note  that:  (i)  the  validity  of  the  last  equality  sign  in  Eq. 
(2.43)  can  be  established  using  Eqs.  (2.2),  (2.3),  (2.5),  (2.7),  and  (2.33);  and  (ii)  because 

dfm/dx-  / dul){dul  / cbc) ,  Eq.  (2.38b)  is  equivalent  to 


->  -» 

du  .du 

- +  A —  =  b 

dt  dx 


(2.44) 


=  / 


(2.45) 


Furthermore,  by  using  Eqs.  (2.23),  (2.33),  (2.35),  (2.36),  (2.42),  and  (2.43),  it  can  be  shown  that 

pu 

2  V  N~l 

f™  +ppp  +  Li= i  PPPi 

puE  +  u{ppp  +  '  p‘Pp, ) 

UP\ 
up2 

uPn-i 

According  to  a  theorem  given  on  pp.l  1-12  in  [22],  Eq.  (2.45)  implies  that  each  flux  function  fm, 
m  =  1,2, ...,N  +  2,  is  a  homogeneous  function  of  degree  1  in  the  variables  w,,w2,...,ww+2  . 


\ 

'  pu  " 

/7M2  +  p 

u(pE  +  p) 

= 

UP\ 

up2 

/ 

v  aPw-i  > 
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For  a  reason  that  will  become  clear  soon,  next  we  will  introduce  a  set  of  new  independent 
variables  um,  m  =  l,2,...,N  +  2,i.e., 


def  def  def  def 

w,  =  p,  U2=U,  u3=e,  W3 +i=Y„  (i  =  l,2,...,N-l) 
By  using  Eqs.  (2.3),  (2.33),  (2.35),  and  (2.46),  it  can  be  shown  that 


(2.46) 


«,=»,,  u2= -  ,  u.s+i=—  =  (2.47) 

Ul  U\  2  M,  )  M, 


«!=«!,  u2=uxu2,  ,  «3+i=“iM3+;  (/  =  1,2,. ..,//-!)  (2.48) 

Let  T  be  the  (N  +2)x(N  +2)  matrix  with  the  element  on  the  mth  row  and  the  f  th  column 
being  3^  / .  Then,  by  the  chain  rule,  it  can  be  shown  that  T"1  (i.e.,  the  inverse  of  T)  is  the 
(N  +2)x(N  +2)  matrix  with  the  element  on  the  mth  row  and  the  l\h  column  being  du^  tdu, . 
By  using  Eqs.  (2.47)  and  (2.48),  one  has 


1 

0 

0 

0 

0  ... 

0 

u 

1 

0 

0 

0  ... 

0 

p 

P 

u2  -E 

u 

1 

0 

0  ... 

0 

P 

P 

p 

_ 

1 

0 

0 

1 

0  ... 

0 

P 

p 

_£± 

0 

0 

0 

1 

0 

p 

. 

• 

. 

p 

. 

Pn-i 

0 

0 

0 

0  ... 

1 
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(f’  =  l>2,...,iV— 1) 


(2.53a) 


dV  dV  a ), 

— '-  +  u — 
dt  d x  p 

From  the  above  discussions,  one  concludes  that  the  system  of  equations  represented  by  Eq. 
(2.38a)  is  equivalent  to  the  system  of  equations  formed  by  Eqs.  (2.50),  (2.51a),  (2.52),  and 
(2.53a).  With  the  aid  of  Eqs.  (2.46)  and  (2.49),  and 


def 

A  = 


0  0 


0  P  u 
P 

0  0  0 
0  0  0 


u 

0 


0 


0  0 


0 

u 


0  N 

0 

0 

0 


0  0  0  0  0 


(2.56) 


The  latter  system  can  be  expressed  as 

dt  dx 

where  u  is  the  column  matrix  formed  by  ,  m  =  1, 2, ...,  N  +  2 . 


(2.57) 


Next,  by  the  chain  rule,  one  has  3^  /dt  =  ^^2(3m„  / dut){dul  / dx) .  Thus,  by  the  definition  of 
T,  one  concludes  that  du/dt  =  Tdu/dt  .  Similarly,  one  has  dul  dx  =  Tdu  /  dx  .  Thus,  by 
multiplying  Eq.  (2.57)  from  the  left  withT-1,  one  concludes  that 


dt  dx 


(2.58) 


Because  Eq.  (2.58)  is  equivalent  to  Eq.  (2.44),  one  concludes  that  (A -T~]AT  }dU/dx  =  0 .  The 
last  identity  is  valid  if  and  only  if 

A  =  T~lAT  (2.59) 

In  fact,  with  the  aid  of  Eq.  (2.23),  Eq.  (2.59)  can  also  be  established  directly  using  Eqs.  (2.42), 
(2.49)  and  (2.56). 

According  to  Eq.  (2.59),  A  and  A  are  related  by  a  similarity  transformation.  Thus  they 
have  identical  eigenvalues,  counting  multiplicity  [22].  However,  compared  with  A,  A  has  a 
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simpler  form.  In  particular,  the  (N  +2)  x  (N  +2)  matrix  A  has  the  special  property  that  the  only 
possible  non-zero  element  on  each  of  its  lowest  N-l  rows  is  the  diagonal  element  u.  As  will  be 
shown  shortly,  because  of  this  special  property,  the  eigenvalues  of  A  (and  therefore  those  oL4) 
can  be  evaluated  in  a  straightforward  fashion. 

To  proceed,  let  M  be  any  KxK  matrix  with  mt  j  being  its  element  at  the  z'th  row 

and  yth  column.  Let  det  M  be  the  determinant  of  M.  Then,  according  to  a  basic  theory  of 
determinant,  for  any  given  i  =  1, 2, ...,  K , 


K 

det  M  =  ^  mi  x  (cofactor  of mi  j )  (2.60) 

j= i 

Furthermore,  let  (i)  IK  denote  the  KxK  identity  matrix  for  any  integer  K  >0;  and  (ii)  Abe  a 
scalar.  Then  the  special  property  of  A  referred  to  in  the  last  paragraph  coupled  with  a  repeated 
application  of  Eq.  (2.60)  implies  that 


det  (A  -  AIn+2  )  =  (u-A)n^  x  det(^3  -  A/3 ) 
where  A^  is  the  3  x  3  submatrix  located  at  the  top-left  comer  of  A ,  i.e., 

f  'v 

up  0 


(2.61) 


a3 


def 


4  «  K 

p  p 

0  P  u 


(2.62) 


Because  the  eigenvalues  of  A  (and  thus  those  of  A)  are  the  roots  ofdet(A-AIN+2)  =  0,  by 
combining  (i)  Eq.  (2.61),  (ii) 


det(^3  -AI3)  =  (u-  A) 


(22- A)2  -  — (1  +  — ) 


(2.63) 


and  (iii)  Eq.  (2.26),  one  concludes  that  A  (and  thus  ,4)  has  three  distinct  eigenvalues,  i.e., 
u  +  af,u-af,  and  u,  with  u  being  an  eigenvalue  of  multiplicity  N. 

Even  though  they  are  not  needed  in  the  CE/SE  development,  for  completeness,  the  eigenvectors 
of  the  matrices  A  and  A  will  be  constructed  and  given  explicitly  in  the  remainder  of  this 
subsection. 
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Let  M‘  denote  the  transpose  of  any  matrix  M.  Let  <j>  =  </>N+2 )  be  an  eigenvector 

of  A  with  the  eigenvalue  u  +  af,  i.e.,  it  is  a  nontrivial  solution  to 


1 - 1 

1^ 

1 

'IT' 

+ 

+ 

1 - 1 

■^1 

II 

o 

(2.64) 

According  to  Eq.  (2.56),  Eq.  (2.64)  is  equivalent  to 

-a A  +M=° 

(2.65a) 

P  p  N- 1 

(2.65b) 

-<fi2-af</>3=0 

P 

(2.65c) 

and 

1-4,5,. ..,N  +  2 

(2.65d) 

Because  af>  0  and  p  >  0 ,  Eqs.  (2.65a),  (2.65c),  and  (2.65d)  imply  that 


tf>\  —  <f>2 j  =  <j>2,  and^=0  for £=4, 5,...JV  +  2 


a, 


paf 


(2.66) 


Moreover,  with  the  aid  of  Eq.  (2.26),  it  can  be  shown  that  Eq.  (2.65b)  is  automatically  satisfied 
if  Eq.  (2.66)  is  assumed.  As  such,  Eq.  (2.64)  is  equivalent  to  Eq.  (2.66),  i.e.,  any  nonzero  (j>  that 
satisfies  Eq.  (2.66)  is  an  eigenvector  of  A  with  the  eigenvalue  u  +  af .  Let^2  =1 .  Then  Eq. 
(2.66)  implies  that  the  (N +  2)  x  1  column  matrix 


-d)  def 

<P  = 


f 


Al,— ,0,0,..,0 


A' 


\af  paf 

is  an  eigenvector  of  A  with  the  eigenvalue!/  +  af ,  i.e., 

[^-(m  +  o/)4+2]?1)=  0 

Similarly,  it  can  be  shown  that  the  (A  +  2)x  1  column  matrix 


,l,-^-,0,0,...,0 


-*(2)  def 

<1  =1  “ 

V  Qf  P°f 

is  an  eigenvector  of  A  with  the  eigenvalue// -af,  i.e., 


V 


(2.67) 


(2.68) 


(2.69) 
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(2.70) 


[A-(u-af)IN+2]/2)  =  0 

Next,  let  (f>  =  (^,,^2,...,^v+2)  be  an  eigenvector  of  ^4  with  the  eigenvalue  u,  i.e.,  it  is  a 
nontrivial  solution  to 

[A-uIN+2]t  =  0  (2.71) 

According  to  Eq.  (2.56),  Eq.  (2.71)  is  equivalent  to 

p<t>2  =0  (2.72a) 

3  =  0  (2.72b) 

and 

^2  =  0  (2.72c) 

P 

Because  p  >  0  and  P  >  0 ,  Eq.  (2.72a)  and  (2.72c)  reduce  to 

(2.73) 

Thus  Eq.  (2.71)  represents  two  independent  conditions  (i.e.,  Eqs.  (2.72b)  and  (2.73))  for  N+  2 
variables.  As  such  it  can  have  only  TV  linearly  independent  solutions.  Let 

<t>  ={<f>AAm\-A, \),  m  =  3,4,...,N  +  2  (2.74) 

where  (i) 

$3)  =  ;  and  $3)  =  Sf  for  £  =  2,3,...,N  +  2  (2.75a) 

P 

and  (ii)  for  anym  =  4,5,...,N  +  2, 

2 

^(m)  = - A  -  and  ^(-)  =  s?  for  t  =  2, 3, . . . ,  N  +  2  (2.75b) 

P 

with  (any  m  and  t  )  being  the  Kronecher  delta  symbol.  Then  it  can  be  shown 
that  A'* ,  m  =  3, 4, . . . ,  N  +  2 ,  are  eigenvectors  of  A  with  the  eigenvalue  u,  i.e., 

[A-uIN+2\fm)  =  0,  w  =  3,4,...,A  +  2  (2.76) 

As  will  be  shown,  these  eigenvectors  are  also  linearly  independent. 
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Let  O  be  the  (N  +  2)  x  (A  +  2)  matrix  formed  using  the  column  matrix  ^ ,  m  =  1,2, . . A  +  2 , 
i.e.,  for  any  l ,  m  =  1, 2, . . . ,  A  +  2 ,  the  element  on  the  t  th  row  and  m  th  column  of  <&  is  the  l  th 


element  of  the  column  matrix  .  Then  (i) 


<l>  = 


af 

1 

P 

paf 

0 

0 

0 


paf 

0 

0 


0 


PPe 

P 

0 

1 

0 

0 

0 


ppp, 

p 

0 


1 

0 

0 


P2P 


Pi 


P2P, 


p 

0 


0 

1 

0 


(ii)  with  the  aid  of  Eqs.  (2.4),  (2.20)  and  (2.26),  one  has 


det  d>  =  2 


f  \ 

P  +  Pt 

V  af  J 


>0 


Pn- i 


P 

0 


0 

0 


i.e.,  <6  is  nonsingular;  and  (iii)  Eqs.  (2.68),  (2.70),  and  (2.76)  are  equivalent  to 


(2.77) 


(2.78) 


(2.79) 


where  A  is  the  (N  +  2)x(N  +  2)  diagonal  matrix  defined  by 


ru  +  af 

0 

0 

0  . 

.  0" 

0 

u-af 

0 

0  . 

.  0 

0 

0 

u 

0  . 

.  0 

0 

0 

0 

u  . 

.  0 

<  0 

0 

0 

0  . 

■ 

(2.80) 


Note  that,  because  <E>  is  formed  using  the  eigenvectors ,m  =  1,2,...,  A  +  2,  nonsingulaity  of 
c£>  implies  that  these  eigenvectors  are  linearly  independent. 

For  any  square  matrix  M  and  any  eigenvalue  a  of  M,  let  ES  (M :  a)  denote  the 

eigenspace  of  M  with  the  eigenvalue  a.  Then  because  (i)  (f)^ ,  w  =  1, 2, . . . , vV  +  2 ,  are  linear 
independent,  and  (ii)  the  multiplicities  of  the  eigenvalues  u  +  af ,  u-af,  and  u  are  1,1,  and  A, 
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respectively,  one  concludes  that  (i)  forms  a  basis  of  ES(A\  a +  \  (ii)  forms  a 
basis  of ES(A:a-ufy,  and  (iii)  forms  a  basis  of ES(A:a  +  u).  In  other 

words,  (i)  a  column  matrix  is  an  eigenvector  of  A  with  the  eigenvalue  u  +  af  if  and  only  it  can 
be  expressed  as  with  /?  being  a  nonzero  scalar;  (ii)  a  column  matrix  is  an  eigenvector  of 
A  with  the  eigenvalue  u-af  if  and  only  it  can  be  expressed  as  p(f>{2)  with  (3  being  a  nonzero 
scalar;  and  (iii)  a  column  matrix  is  an  eigenvector  of  A  with  the  eigenvalue  u  if  and  only  it  can 

be  expressed  as  a  linear  combination  of  <j>^m\m  =  3,4,...,N  +  2 ,  where  not  all  combination 
coefficients  vanish. 

Next,  by  using  the  relation  A  =  TAT~ 1  which  follows  from  Eq.  (2.59),  Eq.  (2.79)  implies 


that 

A*?  =  *¥  A 

(2.81) 

where 

def 

-  7”’a> 

(2.82) 

Because  det  (r-1 )  =  pN+i 

>  0  (see  Eq.  (2.49)),  Eqs.  (2.78)  and  (2.82)  implies  that 

det  T  =  det(E~! )  •  det  $  >  0 

(2.83) 

i.e.,  VF  is  nonsingular. 

Let 

_(!)*/ 
V  = 


P-A(u  +  af)A(H  +  uaf)  A  A ..... 
1°/  af  af  af  af 


-(2)  def 
V  = 


— ,—{u-af),—(H -uaf), 

n  n  J  /-§  ^  _ 


af  af 


af  af 


-(3)  def 

V  = 


PPe  PPeU  PPe 
P  ’  P  ’  P 


PlPe  PlPe 
P  ’  P 


v 

Pn- i 

(2.84a) 

af  ) 

Y 

Pn- i 

(2.84b) 

°f  J 

Y 

Pn-\  Pe 

•*5 

(2.84c) 

"■  ,  j 


and,  form  -4,5, ...,N  +  2, 
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p«)d'f(p2p^  p2pPi p2pPm_  e  m/y3 

V7  ~ V  5  9  9 

p  p  p  p 

PPn-\ P  pm-i  ^  3 \/ 


^  om-j 

'  po  I  ? 


•3  PP^P Pm-i  ~  £*m— 3 


P^2  ’ 


(2.84d) 


™=4,5,...,at+2 


Then,  with  the  aid  of  Eqs.  (2.14),  (2.20),  (2.26),  (2.43),  (2.49),  (2.77)  and  (2.82),  it  can  be 
shown  that  the  nonsingular  matrix  is  formed  using  the  column 

matrices  ,m  =  \,2,...,N  +  2  .  As  such,  (i)  these  N  +2  column  matrices  are  linearly 
independent;  and  (ii)  Eq.  (2.81)  is  equivalent  to 

[^-(M  +  fl/)4+2]^(1)  =0  (2.85a) 

[A-(u-af)IN+2y2)  =  0  (2.85b) 

and 

[A-uIN+2\y}{m)  =  0,  m  =  3,4,...,N  +  2  (2.85c) 

Thus  (i)  jy/^  |  forms  a  basis  of  ES  (A :  a  +  uf  ) ;  (ii)  jy/(2^  j  forms  a  basis  of  ES  [A :  a  -  uf  ) ;  and 
(iii)  |y/(3),^,...,y/(JV+2)|forms  a  basis  of  ES  (A  :a  +  u).  In  the  following,  bases  with  simpler 


structure  will  be  constructed. 

As  a  preliminary,  note  that  Eqs.  (2.26)  and  (2.43)  implies  that 

//-£-A«/)2=z--z0+-)=-- 

Pe  P  Pe  P  P  Pe 

Then,  (i)  Eqs.  (2.84a),  (2.84b),  and  (2.84c)  imply  that 


7,0) def  af  -a) 

0  -  —  y  - 

P 


V 


\,u  +  af,H  +  uaf,—,—,...,^LL 

\  P  P  P  J 


-(2)  f  i  rr  A  A_ 

c7  I  A, Qfyii  uaf,  , 


V 


p  p 


(3)  <*/  p  _(3) 

0  =  -t—y  = 

PPe 


( 


\,u,H- 
V  * 


A(°/)2  A  A  Aw-i 


>  >  >•••» 

P  P  P 


and 


(2.86) 


(2.87a) 

(2.87b) 

(2.87c) 
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_  cm— 3  _  cm— 3  _  cm— 3 

pb\  ,  />d2  >  •  •  •?  P&N-1 


(2.87d) 


,(»)*/-(»)  PPP  -0)_  nft  P  PPm_3 


A' 


e  =y/  =  0,0,- 

Pc  P, 


m  =  4,5,...,N  +  2 

Let©  be  the  matrix  formed  using  the  column  matrices 0^ , m  =  1,2,.. ,,N  +  2,  i.e., 

0 
0 


0: 


1 

1 

1 

0 

0 

u  +  af 

« 

l 

V? 

u 

0 

0 

H  +  uaf 

H-uaf 

H 

P2Pa 

p2pft 

Pe 

Pe 

Pe 

El 

El 

El 

P 

0 

P 

P 

P 

r 

El 

El 

El 

0 

p 

P 

P 

P 

r 

Pn- i 

Pn-  i 

Pn- i 

0 

0 

P 

P 

P 

p2p 


ftf-1 


Pe 

0 


Then,  with  the  aid  of  the  fact  that  'Pis  the  matrix  formed  using  ^m\m  =  1,2,. 
(2.87a)-(2.87d)  imply  that 

©  =  ¥0 


where 


def 

Q  = 


a 


0  0 


0 


0 


a 


0 

0  0 


P 


0 

p 
P  Pe 
0 
0 


0  0 

PPA  pp, 

Pe  Pe 

1  0 

0  1 


ft 


0 

0 

ggfla 

Pe 

0 

0 


0  0 

0  0 

0  0  0  0  0  ...  1 

Because  det  (?  =  p(af)  /(p3pe)>0>  Eqs.  (2.89)  and  (2.83)  implies  that 

det  ©  =  det  ¥  ■  det  Q  >  0 


(2.88) 

) 

,N  +  2,  Eqs. 

(2.89) 

(2.90) 

(2.91) 
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i.e.,  ©  is  nonsingular.  In  turn,  this  implies  that  9^ ,  m  =  1, 2, . . . ,  N  +  2 ,  are  linearly 
independent.  Moreover,  because  AQ  =  QA ,  Eqs.  (2.81)  and  (2.89)  implies  that 

A 0  =  ©A  (2.92) 

i.e.,  Eqs.  (2. 85a)-(2.85c)  remain  valid  if,  for  each  +  2,  if/^  is  replaced  by  9^m\ 

As  such,  one  concludes  that  (i)  j  forms  a  basis  of  ESi^A  :a  +  ufy,  (ii)  forms  a  basis 

of  ES  (A :  a  -  uf  ) ;  and  (iii)  j  9^,9^,.  ..,0<iN+2^  j  forms  a  basis  of  ES  (A :  a  +  u) . 

Note  that,  for  the  special  case  N  =  5,  (i)  [^/?  /  -42a f  j  0(1)  ,\^pl  J2af  j  6{2) ,  and  9® , 

respectively,  reduce  to  the  second,  the  third,  and  the  first  column  matrices  contained  in  the 
matrix  S;  and  (ii)  for  each  m  =  4,  5,  6,  7,  9^  reduces  to  the  m  th  column  matrix  contained  in  S. 
Note  that  the  symbols  pi,p2,...,pN_l  and  af  used  here,  respectively,  are  replaced  by 

CVC2,...,CN_X and  a. 

2.4  Evaluation  of  Thermodynamic/Chemical  Parameters 

Let  (i)  h  and  cpl  be  the  specific  enthalpy  and  the  specific  heat  at  constant  pressure 
of  species  i,  respectively;  and  (ii) 

*  def  def  def 

h,  =  Mfr,  s,  =  M'S,,  cpl  =  Mfp,  (2.93) 

Because  Mt  and  st  are  the  molecular  weight  and  the  specific  entropy  of  species  i,  respectively, 
Eq.  (2.93)  implies  that ht,s, ,  and  cpi are  the  molar  specific  enthalpy,  the  molar  specific  entropy, 
and  the  molar  specific  heat  at  constant  pressure  of  species  /,  respectively. 

Because  the  gas  mixture  is  thermally  perfect  (which  is  a  good  approximation  as  long 
as  each  pt is  in  the  order  of  1  atm  or  less),  for  each  species  /,  A,  =hi{T)andcpl  =cpj(T). 
Furthermore,  in  the  current  study,  each  cpi  (F)  is  approximated  using  a  polynomial  of  T,  i.e., 

c  (T)  JL 

^-  =  taJe  (2-94) 

(=o 
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where  the  coefficients  aa  can  be  found  in  NASA  CEA  program  data  base.  By  using  Eq.  (2.94) 
and  the  relation 

W)=  ljP,{T'W'  +hfi  (2.95) 

where  hf  is  the  molar  heat  of  formation  for  species  i  at  the  standard  conditions 
T  =  Tref  =  293. 15 AT  and  pt  =  1  atm ,  it  can  be  shown  that 


h(T)  au  Tt+ 1 

K  i5  1 


(2.96) 


Here 


a  tf^L-V-^-(r  )'+I 
'5  ^f+i(re/) 


hfi 


(2.97) 


can  also  be  found  in  NASA  CEA  Program.  As  such  h.  can  be  evaluated  using  Eq.  (2.96).  In 

/V 

turn,  h:  and  et  can  be  evaluated  using  the  relations  (i)  ht  =  h,  /  M:  and 
(ii )elshl-pllp,  =  h,  -  Rtr . 

Next,  by  using  (i)  Eqs.  (2.1)  and  (2.93),  (ii)i^  =  MtRt,  and  (iii)  the  thermodynamic 
relation  Yds,  =  dh  -dpt  /  pt ,  one  concludes  that 


T  Pi 


(2.98) 


Let  sf°  denotes  the  value  of  s,  at  /?(  atm.  Then  s°  =  s°  (T)  and  Eq.  (2.98)  implies  that 


V(T)= 

By  using  Eqs.  (2.94)  and  (2.99),  it  can  be  shown  that 

^  =  a,(+a„Inr+^7-' 


where 


- - aJnTref  ~  2.  ~HTf^ 


*—>  0 

(=0  1 


(2.99) 


(2.100) 


(2.101) 


can  be  found  in  NASA  CEA  Program.  As  such  s°  can  be  evaluated  using  Eq.  (2.100). 
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Let  h°  (r)andg(°  (T),  respectively,  denote  the  molar  specific  enthalpy  and  the  molar 

specific  Gibbs  free  energy  of  species  i  at  the  temperature  T  and  the  pressure  =  1  atm.  Then,  by 
definition, 

g°  (T)  =  h°  (T)  -  Tsf  (T)  (2. 1 02) 

Because  ht  is  assumed  to  be  independent  of  pressure,  hf(T)  =  h,(T)  in  current  study.  Thus,  with 
the  aid  of  Eqs.  (2.96)  and  (2.100),  Eq.  (2.102)  implies  that 

=  ?iL  +  _  a  )  _  a  inT  -  V  — (T)'  (2.1 03) 

RJ  T  ,0  '°  £(£  +  !) 

To  proceed,  a  set  of  Nr  chemical  reactions  involving  N  chemical  species  is  represented  by  the 
following  reaction  equations: 

fJvrtAi<=>fjy'A  (r  =  \,2,3,...,Nr)  (2.104) 

i=i  i=i 

Here  (i)  Ai  denote  the  /  th  chemical  species,  and  (ii)  v',  and  v" ,  respectively,  denote  the 
stoichiometric  coefficients  of  reactants  and  products  for  species  i  in  the  rth  reaction.  Let 
(i)  Kfr&ndKbr,  respectively,  denote  the  forward  (i.e.,  left  to  right)  and  the  backward  (i.e.,  right 

def 

to  left)  reaction  rate  constants  of  the  r  th  reaction;  (ii)  nt  =  pt  /  M(  =  the  molar  concentration 

(moles  per  unit  volume)  of  species  i;  and  (iii)  («,)^  denote  the  molar  generation  rate  (moles  per 

unit  time  and  per  unit  volume)  of  species  /  contributed  by  the 
r  th  reaction.  Then  the  law  of  mass  action  implies  that 

(2.105) 

V  e=i  (=] 

Here  Kfr  and  Khr  can  be  evaluated  using  the  Arrhenius  forms,  i.e., 

Kf(T)  =  ApB^V[~^\  and  Klr(T)  =  AbrrBlrexp(  |  (2.106) 

\  ■'V'  y  V  > 

with  Afr ,  Abr ,  Bfr ,  Bbr ,  Efr ,  and  Ebr  being  given  constant  coefficients.  The  values  of 

these  coefficients  for  some  chemical  reactions  are  tabulated  in  Appendices  I  and  II. 

If  the  kinetic  data  of  the  backward  reactions  were  not  available,  Kbr  may  be  evaluated 
using  the  relation 
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(2.107) 


Kbr=Kfr/Keqr 

where  Keqr ,  the  equilibrium  constant  for  the  r  th  reaction,  is  a  function  of  T  and  can  be 
evaluated  using  the  relation 


KeqriT)  = 


(  j 


exp 


/  ;  V  n  r  riJ  rj  rp 
»=1 


(2.108) 


Next,  with  the  aid  of  Eq.  (2.105)  and  the  relations:  (i) 


=M,Yd(nl)r 


(2.109) 


(ii)  «,  =piIMi,  and  (iii)  Kfr  =  Kfr  (T )  and  Kbr  =  Kbr  (T)  ,  r  =  1, 2, . . . ,  Nr ,  each  d>,  can  be 
explicitly  expressed  as  a  function  of  pvp2,...,pN  and  T,  i.e.. 


N, 


i=i 


■  Y" 


*=1 


(2.110) 


As  will  be  shown  shortly,  alternatively  each  a,  can  also  be  considered  as  an  implicit  function  of 
the  conservative  variables  um  ,  m  =  1,2,...,  A +  2  .  In  fact,  evaluation  of  da,  / dum  for  all 
i  =  l,2,...,vV-l  and  m  =  4, 5,..., N  +  2  is  required  for  a  later  application.  As  such,  in  the 
remainder  of  this  section  these  derivatives  will  be  expressed  explicitly  in  terms  of  a  set  of  flow 
and  thermodynamic  variables. 

To  proceed,  first  Eq.  (2.1 10)  is  used  to  obtain 


da. 


dPt 


L  =  Mjy£(v'ri-vri) 


r= 1 


KJrVrk 


V’1  v  w"  N 
A.  V.  t — r 


Pk  (= 


nw  -^nw 


Pk  7=i 


,  i,k  =  \, 2,...,N  (2.111) 


and 


dm. 


dT 


T  =  M^(v'n-vn) 


r=\ 


Kf  n  w'"  ii(".) 

i=\  e=\ 


,  i  =  l,2,..., JV  (2.112) 


with 


.  def  dKr(T)  .  def  dK.  (T) 
K,  =  — f—  and  K.  =  br  K  } 

jr  jrji  br 


dT 


dT 


(2.113) 
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Here  (i)  each  dcol  /  dpk  is  evaluated  assuming  that,  except  pk ,  all  other  independent  variables  in 
the  set  [pvp2,...,pN,T]arQ  held  constant;  and  (ii)  each  dmJdT is  evaluated  assuming  that 
pvp2,...,pN  are  held  constant. 

Next,  because  (i)  p:  =  w3+J ,  i  -  1, 2, ...,  N  - 1 ,  and  (ii)  pN=  p-^  *  1  p,=ux-  1  u( ,  one  has 

dp, 


du„ 


=  C',  1  =  1,2, m  =  l,2,...,N  +  2 


(2.114) 


and 


dPn  = 
du_ 


1  if  w  =  1 
0  if  m  =  2,3 
-1  ifm  =  4,5,. ..,N  +  2 


(2.115) 


Here,  as  in  earlier  cases,  any  partial  derivative  with  respect  to  any  conservative  variable  um  is 
evaluated  assuming  that,  except um,  all  other  variables  in  the  set  {ux,u2,...,uN+2}  are  held 
constant. 

Next,  with  the  aid  of  Eq.  (2.35),  the  first  equation  in  Eq.  (2.40)  implies  that 

de  =— [(w2  -  E)dux  - udu2  +  du3  J  (2.116) 


Also,  with  the  aid  of  Eqs.  (2.13)  and  (2.33),  Eqs.  (2.16)  and  (2.1 16)  imply  that 


dT  = 


PP 


V  2 

\ 

w 

LI 2 

~eN 

y 

N+ 2 


dux  -udu2  +  du3  -  ^  (em_3  ~eN)dun 


m- 4 


(2.117) 


Eqs.  (2.1 17)  implies  that  fis  an  implicit  function  of  ux,u2,...,uN+2  with 


dT  1 


f  -.2 


dux  pcv 


■~eK 


dT 

du~ 


u  ,  dT  1 
- -  and 

PC, 


du3  pcv 


and 


9T  ~  6n  g^-,  m  =  A,5,...,N  +  2 


(2.118) 


(2.119) 


dUn,  P^ 

Note  that:  (i)  each  a>i  is  a  function  of  the  independent  variables  px ,  p2 , . . . ,  pN ,  and  T,  and  (ii)  the 
independent  variables  themselves  are  functions  of  um,  m  =  1,2,...,  A +  2.  Thus  each  coi  can  also 
be  considered  as  a  function  of um,  m-  1,2,..., N  +  2.  Moreover,  by  using  the  chain  rule,  each 


29 


debjldum  (/  =  mdm-  \,2,...,N  +  2)  can  be  obtained  with  the  aid  of  Eqs.  (2.111), 

(2.112),  (2.118),  and  (2.119).  In  particular,  for  /  =  1,2, ...,  — 1  and  m=  4,5,...,N+2,  one 
has 


8cbt  _  yl  deb,  dpk  86),  dpN  dcbi  8T 
*=i  8pk  dum  +  dpN  dum  8T  dum 


deb. 


da b  deb, 

-  + 


dpm- 3  9Pn  dT 


m-3 

N. 


eN  em- 3 


-Ki) 


r~l 


\  PC, 


\ 


FT  (n£y, 
Pn  )  i=\ 


\  ym- 3 


-K, 


br 


Pn  ) 


\  Pm-2 


/=! 


+  M, 


re  -e  . 

"  lE(>'»  ->») 


V 


Pc,  J 


1=1  1=1 


/  =  1,2,...,JV-1;  m  =  4,5,...,iV  +  2 


(2.120) 


2.5  Flow  Equations  in  Multi-Dimensional  Space 


In  the  past,  the  fundamental  behavior  of  cavity  flows  was  known  as  two-dimensional  [8]. 
Equation  (2.121)  shows  the  vector  form  of  the  two-dimensional  flow  equations  in  Cartesian 
coordinates,  including  the  continuity  equation,  the  Navier-Stokes  equations,  the  energy 
equation,  and  one  species  equation: 


SU  SF  dG 

—  +  —  + - 

8t  dx  dy 


8FV  8G 


v  _ 


dx  dy 


=  0 


(2.121) 


where  the  flow  variable  vector 


f 


U  = 


\ 


P 
pu 
pv 
pe 

\PYf) 


the  inviscid  flux  vectors  are 


30 


equation  of  state  of  an  ideal  gas,  p  =  (y  -  l)p  s,  where  y  =  Cp/Cv  is  the  specific  heat  ratio  and  it  is 
a  constant  due  to  the  polytropic  gas  assumption.  In  the  viscous  vectors,  txx,  xXy,  and  Xyy  are 
stress  components,  and  qx  and  qy  are  the  heat  conduction  fluxes  in  x  and  y  directions, 
respectively.  Yf  is  the  mass  fraction  of  fuel.  The  diffusion  velocity  components,  u  and  v  are 
calculated  by  Fick’s  law: 
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(2.122a) 


SYf 

Y fiif  =  —D — — 
f  f  dx 

dYf 

Y  v  f=-D^~ 
f  f  dy 


(2.122b) 


where  D  is  the  mass  diffusivity  of  fuel  in  the  gas  mixture.  The  molecular  viscosity  jo.  is 
calculated  using  Sutherland’s  law  [12]  and  the  Lewis  number  Le  =1  is  assumed  to  calculate  the 
mass  diffusivity  D.  In  numerical  calculations,  the  above  governing  equations  are 
nondimensionalized  by  the  free  stream  conditions,  i.e.,  velocity  components  by  u„,  density  by 
poo,  pressure  by  pooUoo2,  and  the  total  energy  by  pooUoo2-  The  subscript  co  denotes  the  free  stream 
condition.  The  cavity  depth  d  is  used  as  the  length  scale,  and  the  time  scale  is  d/Uoo. 
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3.  The  CESE  Method 


The  CESE  method  is  a  novel  numerical  framework  for  high-fidelity  solution  of  hyperbolic 
conservation  laws.  Originally  developed  by  Chang  and  coworkers  [13-18],  the  tenet  of  the 
CESE  method  is  a  unified  treatment  of  space  and  time  in  calculating  flux  balance.  Contrast  to 
modem  upwind  schemes,  no  Riemann  solver  and/or  a  reconstruction  procedure  is  used  as  the 
building  block  of  the  CESE  method.  As  a  result,  the  logic  and  computational  counts  of  the 
CESE  method  are  simpler  and  more  efficient.  Based  on  the  CESE  method,  computer  programs 
for  solving  unsteady  flows  in  one,  two,  and  three  spatial  dimensions  for  structured,  and 
unstructured  meshes,  and  for  meshes  composed  of  mixed  elements  have  been  developed. 
These  solvers  have  been  parallelized  based  on  domain  decomposition  in  conjunction  of  the  use 
of  MPI.  Since  no  Riemann  solver  is  used,  we  have  straightforwardly  extended  the  CESE 
method  for  flows  with  complex  physical  processes,  including  detonation,  cavitations,  and 
MHD. 

Previously,  various  flow  phenomena  have  been  calculated  by  using  the  CESE  method. 
In  particular,  the  CESE  solver  is  capable  of  calculating  high-speed  compressible  flow  as  well  as 
flows  at  very  low  Mach  numbers  without  applying  preconditioning  to  the  governing  equations. 
The  CESE  method  is  indeed  an  all  speed  solver.  Moreover,  the  CESE  method  is  capable  of 
simultaneously  capture  strong  shock  waves  and  the  acoustic  waves  in  the  same  computational 
domain,  while  the  amplitude  of  the  pressure  jump  across  the  shock  wave  would  be  several 
orders  of  magnitude  higher  than  that  of  the  acoustic  waves. 

The  CESE  method  employed  in  the  present  paper  is  based  on  the  use  of  quadrilateral 
cells  on  the  x-y  plane  [17],  which  was  extended  from  the  original  CESE  method.  Note  that  the 
original  CESE  method  for  two-dimensional  flows  was  designed  based  on  the  use  of  triangular 
cells.  In  the  present  paper,  a  brief  discussion  of  this  particular  extension  of  the  CESE  method  is 
provided.  The  discussions  here  will  be  focused  on  the  space-time  geometry  of  the  CESE 
method.  We  remark  that  the  basic  structure  of  the  CESE  method  can  always  be  easily  grasped 
by  visualizing  the  space-time  geometry  of  conservation  element  (CE),  solution  element  (SE), 
and  how  they  facilitate  the  space-time  integration.  The  detailed  algebraic  equations  of  the 
method,  perhaps,  only  reaffirm  the  structure  of  the  method.  We  of  course  refer  the  interested 
readers  to  the  cited  references  for  all  details. 
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To  proceed,  let  E3  denote  a  three-dimensional  Euclidean  space,  in  which  xj  =  x,  X2=  y, 
and  x3=  t.  Let  V-  be  the  divergence  operator  in  E3,  and  hm  def[(/-  fv)m,(g - gv)m,u J  for 

wi=l,  2,  3,  4,  and  5.  Here  ( ~ gv)m> and  um  are  the  111th  components  of  F-Fv,  G-Gv, 

and  U,  respectively,  in  Eq.  (1).  Aided  by  the  above  definition,  for  each  m=l,  2,...,  5,  the  flow 
equations,  Eq.  (1),  becomes 

Vhm=0,  (3.1a) 

Apply  Gauss’s  divergence  theorem  to  Eq.  (3a)  and  we  get 

•<*  =  <>,  (3.1b) 

As  depicted  in  Fig.  1(a),  S(V)  is  the  boundary  of  the  space-time  region  V  in  E3  and  ds  is  a 
surface  element  vector  pointing  outward.  Equation  (3b)  states  that  the  total  space-time  flux  hm 
leaving  volume  V  through  S(V)  vanishes.  All  mathematical  operations  can  be  carried  out  as 
though  E3  were  an  ordinary  three-dimensional  Euclidean  space.  The  CESE  method  is  designed 
to  accurately  integrate  Eq.  (3b)  to  provide  high  fidelity  results  of  evolving  um  in  the  space-time 
domain. 

The  CESE  method  is  a  family  of  numerical  schemes,  with  the  a  scheme  [13]  as  its 
backbone.  Contrast  to  conventional  finite  volume  methods,  the  CESE  method  has  separate 
definitions  of  conservation  element  (CE)  and  solution  element  (SE)  in  constructing  the 
discretized  equations  for  integrating  Eq.  (3b)  in  the  space-time  domain.  CEs  are 
non-overlapping  space-time  sub-domains  such  that  (i)  the  whole  computational  domain  can  be 
filled  by  the  CEs;  (ii)  flux  conservation  can  be  enforced  over  each  CE  and/or  over  a  union  of 
several  neighboring  CEs;  and  (iii)  inside  a  CE,  flow  discontinuity  is  allowed.  On  the  other  hand, 
SE  are  non-overlapping  space-time  sub-domains  such  that  (i)  SE  do  not  generally  coincide  with 
a  CE;  (ii)  the  union  of  all  SEs  does  not  have  to  fill  the  whole  computational  domain;  (iii)  flow 
variables  and  fluxes  are  discontinuous  across  interfaces  of  neighboring  SEs;  and  (ii)  within  a 
SE,  flow  variable  and  flux  function  are  assumed  continuous  and  they  are  approximated  by  using 
a  prescribed  smooth  function.  In  the  present  paper,  a  first-order  Taylor  series  expansion  is  used. 
Thus  the  discretized  flow  variables  and  fluxes  are  linearly  distributed  inside  each  SE. 

The  time  marching  strategy  in  the  CESE  method  is  designed  based  on  a  space-time  staggered 
mesh  stencil  composed  of  CEs  and  SEs.  Note  that  when  integrating  Eq.  (3b)  over  the  boundary 
of  a  CE,  the  surface  element  of  S(V)  in  Eq'.  (3b)  always  lying  inside  a  SE,  where  flow  variables 
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and  fluxes  are  continuous.  We  remark  that  the  paradigm  of  the  Godunov  schemes  is  that  one  has 
to  resort  to  the  use  of  a  Riemann  solver  to  reckon  the  nonlinear  flux  function  at  the  cell 
interfaces.  In  the  CESE  method,  however,  flow  information  propagates  only  in  one  direction 
across  all  cell  interfaces,  i.e.,  towards  future.  Thus,  the  space-time  flux  integration  can  be 
straightforwardly  carried  out  without  reconciling  the  values  of  flux  functions  at  cell  interfaces 
through  the  use  of  a  Riemann  solver.  In  other  words,  contrast  to  the  upwind  methods,  there  is 
no  cell  interface,  across  which  two  way  traffics  of  flow  information  propagate.  Thus,  the  CESE 
method  can  capture  shocks  without  using  a  Riemann  solver.  In  what  follows,  we  discuss 
specific  space-time  geometry  of  the  CE  and  SE  in  integrating  Eq.  (3  b). 

Consider  Fig.  1(b).  The  x-y  plane  is  divided  into  non-overlapping  quadrilaterals.  Two 
neighboring  quadrilaterals  share  a  common  side.  Vertices  and  centroids  of  quadrilaterals  are 
marked  by  dots  and  circles,  respectively.  Q  is  the  centroid  of  the  quadrilateral  B]B2B3B4.  A],  A2, 
A3,  and  A4,  respectively,  are  the  centroids  of  the  four  quadrilaterals  neighboring  to  the 
quadrilateral  B1B2B3B4.  Q*,  marked  by  a  cross  in  Fig.  1(b),  is  the  centroid  of  the  polygon 
AiBiA2B2  A3B3A4B4.  Hereafter,  point  Q*,  which  generally  does  not  coincide  with  point  Q,  is 
referred  to  as  the  solution  point  associated  with  Q.  Note  that  points  A*,A2,A‘, and  A’ ,  which 
are  also  marked  by  crosses,  are  the  solution  points  associated  with  the  centroids  Aj,  A2,  A3,  and 
A4,  respectively. 

To  proceed,  we  consider  Fig.  1(c).  Here  t=  nAt  at  the  nth  time  level,  where  n  =  0, 1/2,  1, 
3/2, . . .  For  a  given  n>0,  Q,  Q’,  and  Q”,  respectively,  denote  the  points  on  the  nlh,  the  (1-1/2  n)*, 
and  (l+l^n)*  time  levels  with  point  Q  being  their  common  spatial  projection.  Other 
space-time  mesh  points,  such  as  those  depicted  in  Fig.  1(c),  and  also  those  not  depicted,  are 
defined  similarly.  In  particular,  Q* ,  A[,A*2 ,  A\ ,  and  A\  lie  on  the  nth  time  level  and  they  are  the 
space-time  solution  mesh  points  associated  with  points  Q,  A],  A2,  A3,  and  A4. 
Q* ,A* ,A2,A3' ,and  A4  lie  on  the  (l-l^n)*  time  level  and  are  the  space-time  solution  mesh 

points  associated  with  points  Q ’,  Av  A2 ,  A3 ,  and  A4 . 

With  the  above  preliminaries,  we  are  ready  to  discuss  the  geometry  of  the  CE  and  SE 
associated  with  point  Q*,  where  the  numerical  solution  of  the  flow  variables  um  at  nth  time  level 
are  calculated  based  on  the  known  flow  solution  in  all  points  at  a  previous  time  level,  i.e.,  n-1/2, 
denoted  by  superscript  prime.  First,  the  solution  element  of  point  Q*,  denoted  by  SE(Q*),  is 
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defined  as  the  union  of  the  five  plane  segments  QQ'B\B[ ,  Q'Q'B'2B'2  ,  Q’Q'B"3B’3 ,  QQ'B'4B4 , 
and  AlB1A2B2A3B3A4B4  and  their  immediate  neighborhoods. 

To  integrate  Eq.  (3b),  four  basic  conservation  elements  (BCEs)  of  point  Q*  are 
constructed,  and  they  are  denoted  by  BCEi(Q),  with  1  =  1,  2,  3,  and  4.  These  four  BCEs  are 
defined  to  be  the  space-time  cylinders  AlBiQB4AlBiQ'  B4  ,  A2B2QB1A'2B'2Q' B\  , 
A3B3QB2A3B3Q  B2 ,  and  A4B4QB3A4B4Q  B3  ,  respectively.  In  addition,  the  compounded 
conservation  element  of  point  Q,  denoted  by  CE(Q),  is  defined  to  be  the  space-time  cylinder 
AiBiA2B2A3B3A4B4AlBiA2B2A3B3A4B4 ,  i.e.,  the  union  of  the  above  four  BCEs. 

To  proceed,  the  set  of  the  space-time  mesh  points  whose  spatial  projections  are  the 
centroids  of  quadrilaterals  depicted  in  Fig.  1(b)  is  denoted  by  D.  and  the  set  of  the  space-time 
mesh  points  whose  spatial  projections  are  the  solution  points  depicted  in  Fig.  1(b)  is  denoted  by 
£2*.  Note  that  the  BCEs  and  the  compounded  CEs  of  any  mesh  point  eO  and  the  SE  of  any 
mesh  point  efi*  are  defined  in  a  manner  identical  to  that  described  earlier  for  point  Q  and  Q*. 
With  the  clear  definitions  of  the  CE  and  SE  in  above,  the  numerical  integration  of  the 
space-time  flux  balance,  i.e.,  Eq.  (3b),  in  the  present  modified  CESE  method  can  be 
summarized  as  follows. 

For  any  Q*  e  fi*  and  ( x,y,t)s  SE(f2*),  the  flow  variables  and  flux  vectors,  i.e., 
uAx>y, 0,  fm{x,y,t) ,  and  gm(x, y,t),  are  approximated  to  their  numerical  counterparts,  i.e., 
um(x,y,t),  and  g'm(x,  yj),  by  using  the  first  order  Taylor  series  expansion  with 

respect  to  Q*  (xQ.,yQ.,tn) .  Thus  the  space-time  flux  vector  h m(x,y,t),  can  be  replaced  by 

h  'm(x,y,t;Q')  and  the  numerical  analogue  of  Eq.  (3b)  for  each  m=l,2,...,5,  is 

<{  fa  *  *  c/s =  0,  (3  2) 

Equation  (3.2)  states  that  the  discretized  total  flux  of  h^,  leaving  CE(Q)  through  its  boundary 

vanishes.  We  note  that  Eq.  (4)  can  be  written  in  terms  of  independent  discrete  solution  variables 
which  are  the  main  flow  variables  and  their  spatial  derivatives,  ( u •mx)Q.  and  (umy)Q.  in  this 

approximation  process. 


36 


By  conducting  the  integration  of  Eq.  (4)  over  the  CE=  BCE1+BCE2+BCE3+BCE4,  the 
discrete  flow  variables  ( um)Q .  associated  with  the  space-time  point  Q*,  can  be 

straightforwardly  evaluated.  This  is  achieved  by  the  aid  of  the  geometrical  information  of 
CE(Q*)  as  shown  in  Fig.  1,  and  the  linear  distribution  of  (um)Q.  in  each  SE  due  to  the  adopted 

first-order  Taylor  series  expansion  for  the  flow  variables  inside  the  SE  as  the  discretization 
process. 

The  calculation  of  the  gradient  variables,  i.e.,  (am)g.  ,  and  (w^)^.  is  base  on  a 

finite-difference  approach  in  conjunction  with  the  standard  CESE  artificial  damping  functions, 
i.e.,  the  a-s-a  scheme,  in  which  parameter  s  is  associated  with  the  overall  damping  effect  and  a 
is  for  shock  capturing.  In  contrast  to  the  original  CESE  method,  the  calculation  of  ( am)  . ,  and 

( umy)Q •  has  nothing  to  do  with  the  space-time  flux  conservation. 

To  calculate  unsteady  flows,  the  non-reflecting  boundary  condition  treatment  is 
critically  important.  Without  an  effective  treatment,  the  reflected  waves  would  inevitably 
contaminate  the  evolving  flow  solutions.  Numerical  treatments  to  achieve  non-reflecting 
boundary  condition  in  the  setting  of  conventional  CFD  methods  have  been  an  active  research 
subject  for  a  long  time.  In  general,  most  of  treatments  were  developed  based  on  theorems  of  the 
partial  differential  equation,  and  they  could  be  categorized  into  the  following  three  groups:  (i) 
applying  the  method  of  characteristics  to  the  discretized  equations,  (ii)  the  use  of  the  buffer 
zone  or  a  perfectly  matched  layer,  and  (iii)  applying  asymptotic  analytical  solution  at  the  far 
field. 

In  the  setting  of  the  CESE  method,  we  only  concern  the  integral  equation.  The  above 
ideas  of  treating  non-reflective  boundary  are  not  applicable.  Instead,  the  non-reflecting 
boundary  condition  treatments  in  the  setting  of  the  CESE  method  is  based  on  flux  conservation 
in  the  vicinity  of  the  computational  boundary  [18].  In  other  words,  the  present  nonreflecting 
boundary  condition  treatment  is  equivalent  to  letting  the  incoming  flux  from  the  interior 
domain  to  the  boundary  CE  smoothly  exit  to  the  exterior  of  the  domain.  In  the  setting  of  the 
CESE  method,  the  numerical  implementation  of  this  flux-based  method  is  extremely  simple 
due  to  the  fact  that  all  flow  information  must  propagate  into  the  future.  Chang  and  coworkers 
[18]  have  provided  detailed  discussions  of  various  implementations  of  the  above  principle.  It 
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has  been  demonstrated  that  only  negligible  reflection  occurs  when  a  shock  passes  through  the 
domain  boundary.  Moreover,  along  a  wall  boundary,  a  unified  boundary  condition  for  viscous 
flows  is  used.  Based  on  local  space-time  flux  conservation,  a  no-slip  condition  will  be 
automatically  enforced  when  the  viscosity  is  not  null.  Again,  the  basic  principle  is  based  on  the 
space-time  flux  conservation  over  CEs  near  the  computational  boundary. 
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4.  Results  and  Discussions 


Three  sets  of  numerical  results  are  presented:  (i)  supersonic  cavity  flow  in  the  supersonic 
combustion  facility  by  Gruber  et  al.  [1, 2]  and  Baurle  et  al.  [5],  (ii)  Stalling  and  Wilcox’s  cavity 
flow  test  [20],  and  (iii)  cavity  flows  with  fuel  injection.  The  first  test  is  to  assess  numerical 
accuracy  of  the  calculated  frequencies.  The  results  will  be  compared  with  previously  reported 
data  and  Rossiter’s  empirical  relation  [6].  The  second  test  is  to  assess  numerical  accuracy  of 
time-averaged  amplitudes  of  pressure  fluctuations  along  the  cavity  walls.  The  results  will  be 
compared  with  the  experimental  data  [20].  In  the  third  test,  cavity  flows  with  downstream  as 
well  as  upstream  injections  are  simulated.  We  will  show  that  a  cavity  flow  with  a  downstream 
transverse  injection  can  effectively  generate  strong  vortices  and  acoustic  waves  for  fuel/air 
mixing  enhancement. 

4.1  Frequency  Calculations 

The  first  numerical  example  follows  the  testing  condition  in  the  US  AFRL  supersonic 
combustion  facility  reported  in  [1,  2,  5].  A  supersonic  flow  at  Mach  2  and  Reynolds  number  of 
4x1 05  passes  a  swallow  cavity  with  L/d  =  7.76,  where  L  and  d  are  the  length  and  depth  of  the 
cavity,  respectively.  The  computational  domain  outside  of  the  cavity  is  0  <  x  <  1 1 .52 ,  and 
0  <  y  <  3.82 ,  where  x  and  y  are  nondimensionalized  by  d.  Mesh  points  were  clustered  at  the 
forward  and  aft  bulkheads,  at  the  plane  spanning  over  the  cavity  mouth,  and  along  the  lateral 
sidewalls  of  the  cavity.  143,000  quadrilateral  elements  are  used  for  the  computational  domain. 
Figure  2  shows  the  mesh,  in  which  one  of  every  five  mesh  lines  is  displayed.  The 
non-reflecting  boundary  condition  is  applied  to  the  free  stream  surfaces  and  outlet.  Initially, 
velocities  inside  the  cavity  are  set  to  zero,  and  the  density  and  pressure  of  the  whole  domain  are 
set  to  the  free  stream  values.  The  time  step  was  chosen  such  that  CFL  «  0.8  based  on  the  free 
stream  condition. 

Figure  3  shows  a  series  of  snapshots  of  pressure  contours,  vorticity  contours,  and 
numerical  Schlieren  images,  which  are  contour  plots  of 

N= 
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These  figures  demonstrate  very  complex  flow  features,  including  traveling  acoustic  waves, 
vortex  generation  at  the  leading  edge,  shedding  vortices  in  the  free  shear  layer,  and  pressure 
waves  impinging  on  and  rebounding  from  the  aft  wall.  The  interactions  between  the  rebounding 
pressure  waves  and  shed  vortices  form  a  feed  back  loop  which  leads  to  self-sustained 
oscillations  as  illustrated  by  Rossiter  [6].  In  Fig.  3(b),  periodic  shear  layer  deflections  in  the 
transverse  direction  could  be  clearly  discerned.  Inward  deflection  results  in  mass  addition  into 
the  cavity;  outward  deflection  expels  mass  from  the  cavity.  This  periodic  mass 
addition/expulsion  mechanism  enhances  fuel/air  mixing.  Moreover,  flapping  shock/expansion 
waves  emanating  from  the  upstream  bulkhead  of  the  cavity,  shown  in  Fig.  3(c),  can  also 
enhance  fuel/air  mixing. 

Figure  4  shows  pressure  histories  on  the  aft  wall  and  on  the  floor.  The  flow  has  reached 
a  self-sustained  oscillatory  state  after  about  15  tc,  where  tc  =  d /  U*,.  However,  the  oscillation 
pattern  changes  from  cycle  to  cycle,  and  we  cannot  clearly  identify  the  period  of  the  oscillation 
cycles.  This  is  consistent  with  experimental  observation  reported  in  [19].  The  amplitude  of  the 
pressure  oscillations  at  the  aft  wall  is  much  higher  than  that  at  the  cavity  floor  due  to  the  mass 
addition/expulsion  mechanism  near  the  aft  wall.  Figure  5  shows  the  frequency  spectra  of  the 
pressure  data  in  Fig.  4.  The  predicted  values  of  the  dominant  frequencies  compare  well  with  the 
Rossiter  relation  [6]  and  the  numerical  results  by  Baurle  et  al.  [5],  The  Rossiter  formula  is 


(4.2) 


where  fm  is  the  resonant  frequency  corresponding  to  the  /nth  mode,  U  is  the  main  stream 
velocity,  L  is  the  cavity  length,  =  0.513,  and  K=  0.57  [5].  Figure  5  shows  five  calculated 
dominant  frequencies,  i.e.,  2523,  3533,  4290,  5804,  and  8579Hz,  compare  well  with  the  3rd,  4th, 
5th,  6th,  and  9th  mode  predicted  by  the  Rossiter  relation.  The  first  and  second  frequency  modes 
also  compared  well  with  Baurle’s  calculation. 

To  proceed,  we  perform  simulation  of  the  same  cavity  flow  with  an  added  transverse 
injection  at  upstream  of  the  cavity.  The  free  stream  flow  condition  and  the  cavity  geometry  are 
identical  to  that  shown  in  Figs.  3-5.  The  injection  jet  opening  is  0.2  d,  and  its  center  is  located 
1 .0  d  upstream  from  the  leading  edge  of  the  cavity.  A  choked  jet  with  a  uniform  condition  at  the 
opening  is  imposed: 
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